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Abstract. - The elastic properties of two-component bilayer membranes are studied using a 
coarse grain model for amphiphilic molecules. The two species of amphiphiles considered here 
differ only in their length. Molecular Dynamics simulations are performed in order to analyze 
the shape fluctuations of the two-component bilayer membranes and to determine their bending 
rigidity. Both the bending rigidity and its inverse are found to be nonmonotonic functions of 
the mole fraction xb of the shorter B-amphiphiles and, thus, do not satisfy a simple lever 
rule. The intrinsic area of the bilayer also exhibits a nonmonotonic dependence on xb and a 
maximum close to xb — 1/2. 



Biological membranes are multicomponent systems consisting of mixtures of many different 
lipids and proteins. Because the composition of lipid bilayers affects their physical properties 
and biological functions, this composition varies from one organism to the other and between 
organelles in the same cell [1]. Biological membranes contain a fluid bilayer which is highly 
flexible and, thus, can easily change its shape. Typical examples are provided by the plasma 
membranes of red and white blood cells which are so flexible that they can move through 
rather small capillaries. This flexibility is also responsible for the thermally excited shape 
fluctuations of biomembranes in physiological conditions, the amplitude of which depends on 
temperature, membrane composition and mechanical constraints. One important example is 
provided by mixed membranes containing phospholipids and cholesterol which exhibit a strong 
increase in rigidity with increasing amount of cholesterol [2-4]. In erythrocytes, on the other 
hand, the amplitude of the shape fluctuations is influenced by the spectrin-ankyrin network 
which is coupled to the interior of the plasma membrane [5,6]. The study of these fluctuations 
can thus provide information on the elastic properties of the bilayer membrane and on how 
these properties depend on membrane composition. One important elastic parameter is the 
bending rigidity [7], which describes the resistance of the membrane to bending, and which 
can be obtained from the spectrum of the shape fluctuations [8] . 
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In spite of their complex chemical composition, all biomembranes have the same basic 
structure: a bilayer of lipid and protein molecules. Therefore the simplest model for a bio- 
logical membrane is a bilayer of amphiphilic molecules as studied in the present work. The 
bending rigidity of such bilayers has already been determined via computer simulations both 
for coarse grain models of amphiphile- water systems [9,10] and for models with atomic reso- 
lution [11]. In all of these previous studies, the bilayers contained only a single component in 
contrast to natural membranes which contain many different types of lipids and proteins. 

In the present article, we study bilayer membranes with two components and determine the 
dependence of their elastic properties on membrane composition. We use a coarse grain model 
for the amphiphiles and for the water particles which we investigate by Molecular Dynamics 
simulations. The two species of amphiphiles considered here differ only in their length. This 
choice allows us to focus on the contribution from the mismatch of the amphiphile tails. 

Two types of theories have previously addressed the bending elasticity of two-component 
bilayers. First, a self-consistent (or molecular field) theory was used to estimate the contri- 
bution to the bending rigidity arising from the conformations of the amphiphilic tails [12]. 
Secondly, the so-called hat model for local curvature fluctuations was used in order to de- 
rive an expression for the bending rigidity of two-component bilayers [13]. The predictions 
of these two theories are different as far as the functional dependence of bending rigidity on 
composition is concerned. On the one hand, the results of the self-consistent theory for the 
chain conformations depend on the types of packing constraint acting on the conformations 
of a single tail. The hat model, on the other hand, predicts that the inverse bending rigidity 
(or flexibility) satisfies a simple 'lever rule' corresponding to a linear interpolation between 
the two limiting values for the one-component bilayers. 

In contrast to the analytical but approximate theories in [12] and [13], our simulations 
take all types of fluctuations into account: different tail conformations; curvature (or bending) 
fluctuations; molecular protrusions; and composition fluctuations arising from lateral diffusion. 
We find that both the bending rigidity and its inverse exhibit a nonmonotonous behavior as 
a function of the mole fraction xb ■ Our results are qualitatively similar to the self-consistent 
theory, provided one uses the packing constraint of 'constant area' per molecule, but show 
that a simple lever rule for the inverse bending rigidity does not apply in general. 

The coarse grain model used here consists of several types of particles: water particles, 
which correspond to several water molecules, as well as head group and tail particles, which 
represent the head groups and the hydrocarbon tails of the amphiphilic molecules, respectively. 
The two species of amphiphiles will be denoted by A and B. Both types of amphiphiles are 
composed of one head group particle, which is hydrophilic, and of several tail particles, which 
are hydrophobic. The different length of the two molecular species arises from the different 
number of tail particles which is four and two for A and B, respectively, see figure ^ The 




Figure 1 - Cartoon of a two-component bilayer. 



particles experience attractive and repulsive pair-interactions which are modelled by Lennard- 
Jones (LJ) and Soft-Core potential (SC), respectively. The parametrization of these potentials 
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has been previously described in [14] where the same model was used in order to study the 
lateral and transverse diffusion within the two-component bilayer membranes. Within this 
model, all quantities can be expressed in terms of three basic scales: the particles mass to, 
the LJ interaction energy e and the LJ radius p. As in [14] and [15], we chose the numerical 
values vtlNav — 3.6 • 10~^kg, eiV^i^ = 2kJ, and p = 1/3 nm, where Nav is the Avogadro 
number. Our MD simulations are carried out for a cuboidal box of constant volume and 
periodic boundary conditions using a constant temperature algorithm with kgT = 1.35e. The 
starting configurations for the two-component bilayer simulations are such that each monolayer 
contains the same number of A and B molecules. We checked that, in each monolayer, the 
average number of B molecules did not change during the simulations. While multicomponent 
bilayer membranes may exhibit deviations from ideal mixing behavior [16] and may, in general, 
undergo phase separation, the two-component bilayers studied here stayed in the one-phase 
region over the whole range < xb ^ 1 of mole fractions for the B molecules. Thus, our 
bilayers correspond to a lipid mixture above the liquidus-solidus line. 

The elastic properties of a fluid membrane which has vanishing spontaneous curvature 
are governed by its bending rigidity k and by its surface tension a. These quantities can 
be determined from the spectrum of the shape fluctuations. In order to do so, we choose 
the (a;, 2/)-plane to be parallel to the bilayer membrane. The shape of the membrane is then 
described by the height function h{x^ y) which measures the distance of its midsurface from the 
{x, y)-plane. We define the Fourier coefhcients = A'^^ dxdy exp [—i{xqx + y^ly)] h{x, y), 
with q = [qx, Qy), where Ap is the projected area of the membrane. Its fluctuation spectrum, 
S{q) = ( |/iq|^), depends only on q = |q| and exhibits the functional form 



for short-wavelength molecular protrusions [17]. Apart from the temperature T and the pro- 
jected area Ap, the two spectra as given by and ((SJ contain three parameters: the surface 
tension ct, the bending rigidity k, and the protrusion tension (jpr- As discussed in [15], the 
lateral size of the simulation box can be adjusted in order to obtain a bilayer with vanish- 
ing thermodynamic tension. This latter procedure consists in determining, for fixed particle 
number, fixed volume, and fixed temperature, the stress (or pressure) tensor which has a tan- 
gential component, St, and a normal component, E^r [15]. Both components depend on the 
coordinate z perpendicular to the membrane. The thermodynamic tension S is then obtained 
from Ti = J dz[J^T{z) — "Sn^z)]. It is intuitively appealing to identify the thermodynamic 
tension S with the surface tension a which governs the fluctuation spectrum. Indeed, this 
assumption was implicit in our previous work in [9] and was found to be satisfied within the 
accuracy of the simulation data obtained there. However, more extensive simulations have 
shown that this identity does not hold in general [18]. This difference is related to the fact that 
E is conjugate to the projected area, whereas a is conjugate to the intrinsic area, a distinction 
which has been previously discussed, e.g., in [19]. Therefore, in the present work, we used an 
improved procedure in order to determine the state with a — 0: for fixed particle number and 
fixed temperature, we adjusted the lateral box size in order to obtain a fluctuation spectrum 
S{q) which can be well fitted, for small q, by the form iQJ with ct = [18]. 

First, let us discuss this fitting procedure for one-component bilayers corresponding to 
vanishing mole fraction xb = 0. An example of the measured fluctuation spectrum S{q) for 
TVa = 1152. is shown in figure |21 The lines, which represent the functional forms as given by 




(1) 



for long- wavelength bending modes [8] , and the somewhat different form 




(2) 
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with (7 — and by respectively, are plotted in the same figure. From these fits, we 
obtain the values k = 3.0 ± O.iksT and = 2.50 ± 0.03 e/p^. Essentially the same values 
are found for Na = 512: k = 3.20 ± 0.06 /fcsT and CTp,. = 2.65 ± 0.05 e/p^ [18]. For surface 




Figure 2 - Fluctuation spectrum S for one and two-component bilayers with A'^a = 1152 (a;B = 0) 
and A*' = 512 (xb = 0.4), respectively, in tension-free state. The lines are obtained by fitting the 
fluctuation spectra of the tensionless systems, using eq. Q (full lines) and eq. Q (dashed lines) at 
small and large q, respectively. The units are given in brackets. 

tension cr = 0, the intrinsic area per molecule, a.m, is expected to attain a certain value which 
is determined by the optimal packing of the molecules within the bilayer and, thus, should 
not depend on the overall size of the membrane. Indeed, if we adjust the box size to obtain 
vanishing surface tension cr ~ 0, our MD simulations lead to ai„ = 2.26 ±0.01 for Np^ — 512 
and to ain — 2.277±0.005 for Np^ — 1152 which are identical within the numerical accuracy. 

Next, we consider two-component bilayers with mole fraction < xb < 1. To save 
computation time, we varied the mole fraction for fixed total number N — TVa + A'b of the 
amphiphiles with N = 512. In order to attain tensionless states of the two-component bilayers, 
the projected area Ap has to be gradually decreased as xb is increased, since the projected 
area per molecule, ap = 2Ap/N , is smaller for the B-amphiphile than for the A-amphiphile. 
The measured fluctuation spectrum S{q) for xb — 0.4 is also shown in figure |2] 

In figure |3fa), the bending rigidity k is plotted as a function of a;B- For the pure A-bilayer 
and the pure B-bilayer, we find ka — 3.20 ± 0.06 fcsT and kb = 1.3 ± 0.1 fc^T, respectively. 
As shown in figure |2Ia), the bending rigidity is nonmonotonic for intermediate values of xb 
and exhibits a minimum for x^ ~ 0.6 with k — 1.0 ±0.1 fcsT. Figure|3Jb) displays the inverse 
bending rigidity l//t as a function of the mole fraction xb • For this latter quantity, a simple 
lever rule as given by 1/ k = x^/ + [1 — x^)/ k\ has been proposed in [13]. Inspection of 
figure|3Ib) shows that such a lever rule does not apply here. However, our data are consistent 
with the more general proposal that the bending rigidity behaves smoothly both at xb = 
and at Xb = 1. Indeed, we conclude from the data in figure |3Jb) that 

l/«;« (l/«;A)(l + </'ABa;s) with <?!)ab = 3.5 ± 0.1 (3) 

for small xb and that 

1/k« (l/KB)(l + 0BAa;A) with 0BA = 1.O±O.2 (4) 



for small xa = 1 — a^B . 
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Figure 3 - (a) Bending rigidity k and (b) inverse bending rigidity of two-component A- B-bilayers vs. 
the shorter amphiphile mole concentration xb ■ The error bars represent statistical errors obtained by 
the fit of the fluctuation spectrum data to eq. Q. In (a) the line is obtained by fitting the data with 
a smoothing spline. In (b) the full and dotted line correspond to linear fit to the data for xb < 0.4 
and Xb > 0.8, respectively. The units are given in brackets. 



We also measured the intrinsic bilayer area A, defined hy A = dxdy^yi + (V/i)^. In 
Fig. 21 the intrinsic area per amphiphile, aj„ — 2A/N , is plotted as a function of x-b for ten- 
sionless bilayer states. As shown in this figure, the intrinsic molecular area ai„ is again found 
to be a nonmonotonic function of x^: for x^ — 0, it has the value ~ 2.26p^, but exhibits 
a maximum at x-q — 0.4 with ~ 2.39p^. The change in the intrinsic area is obviously con- 

2.42 n ~ ~ ~ ~ n 
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Figure 4 - Area per amphiphile ain ~ 2A/N vs. the shorter amphiphile mole concentration xb- The 
error bars represent standard errors obtained from the area fluctuations. The line is a guide to the 
eye. The units are given in brackets. 

nected to the change in the bending rigidity. The decrease of the latter quantity corresponds 
to a more flexible bilayer with larger shape fluctuations and a corresponding increase in the 
intrinsic area. Inspection of flgure 0] shows that the pure A-bilayer has essentially the same 
intrinsic molecular area as the pure B-bilayer, i.e., ain{x-B = 0) ~ ain{x-B = 1) ~ 2.26/3^. In 
contrast, the corresponding projected molecular areas Op are found to be ap{xB = 0) = 2.10p^ 
and ap{xB = 1) = 1.96p^ which reflects the increased shape fluctuations for xb — 1- 
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Finally, the protrusion tension apr as introduced in ^ also exhibits a nonmonotonic 
behavior as shown in Fig. [S] A decrease in the protrusion tension implies an increase in the 
amplitude of the local protrusions which become energetically more favorable. As shown in 
Fig. apr decreases with increasing xb up to xb = 0.2, then stays essentially constant up to 
Xb = 0.8, and finally increases again up to xb = 1. This indicates that a small amount of 
shorter B-amphiphiles is sufficient to roughen the bilayer surface, and further addition of these 
amphiphiles has essentially no effect on the local protrusions until one reaches another regime 
characterized by a bilayer of B-amphiphiles with a few longer A-amphiphiles. Comparing 
the results for bending rigidity, intrinsic area and protrusion tension, one concludes that 
the changes in the intrinsic area depend both on the shape fluctuations and on the bilayer 
roughness arising from the local mismatch between the two types of amphiphiles. 

The bending rigidity of amphiphilic bilayers has been measured in many experiments: while 
phospholipid membranes [4,8,20-22] are characterized by a bending rigidity of the order of 
tens of ksT, lamellar and fluid microemulsion phases composed of single chain surfactants, 
short chain alcohols, oil and water [23-26] exhibit a bending rigidity of the order of fc^T. 
The bilayer studied in our Molecular Dynamics simulations are composed of single-chain 
amphiphiles and have bending rigidities which are of the same order of magnitude as those 
studied in the second group of experiments. In Refs. [23-25], the response of lamellar systems 
to the insertion of shorter cosurfactants was studied experimentally. In these systems, the 
main surfactant was sodium dodecyl sulfate (SDS), which is a single chain surfactant, while 
the cosurfactants were alcohol molecules: it was found that the bending rigidity decreases as 
the fraction of shorter chain cosurfactant increases. In Ref. [25], short surfactant (pentanol) 
molecules were added to bilayers composed of two-chain amphiphiles (DMPC) and a decrease 
of the bending rigidity was observed. Two-component bilayer membranes were also studied 




Figure 5 - Protrusion tension apr vs. the shorter amphiphile mole concentration xb- The error bars 
represent statistical errors obtained by the fit of the ffuctuation spectrum data to eq. ||5J. The fine 
is a guide to the eye. The units are given in brackets. 



theoretically by Szleifer et al. [12] using a simple mean field theory to treat chain packing. 
In the framework of these theories, which are based on smooth bending configurations and, 
thus, neglect protrusions, the longer amphiphiles obtain more configurational freedom as the 
shorter amphiphiles are added which leads to an increased flexibility of the bilayers. The 
behavior of the bending rigidity as found in our simulations for small xb is consistent with 
such a mechanism. However, we also find that both the bending rigidity and the protrusion 
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tension decrease as one adds long A-amphiphiles to a bilayer of short B-amphiphiles for large 
values of xb, see Figs. |21and[Sl In this latter regime, the addition of the A-amphiphiles 
induces more protrusions which act to reduce the bending rigidity [17]. 

In conclusion, we have presented the results of Molecular Dynamics simulations for a coarse 
grain model of two-component bilayers. We show that the functional dependence of the in- 
verse bending rigidity on the membrane composition does not follow a simple lever rule. In 
contrast, we find that both the addition of short amphiphiles to a bilayer of long ones and 
the addition of long amphiphiles to a bilayer of short ones leads to an increase in molecu- 
lar protrusions and to a decrease in the bending rigidity. We also observe a nonmonotonic 
functional dependence of the intrinsic area on the membrane composition which is strongly 
correlated with the behavior of both local protrusions and long-wavelength shape fluctua- 
tions. The simulation approach described here can be extended to bilayer membranes with 
three components which have recently been shown to lead to coexisting membrane domains or 
'rafts' and to more realistic models of biological membranes which contain membrane proteins. 
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